GWAS combined with QTL mapping reveals the genetic loci of leaf morphological characters in Nicotiana tabacum

Background Leaf morphology plays a crucial role in photosynthetic efficiency and yield potential in crops. Cigar tobacco plants, which are derived from common tobacco (Nicotiana tabacum L.), possess special leaf characteristics including thin and delicate leaves with few visible veins, making it a good system for studying the genetic basis of leaf morphological characters. Results In this study, GWAS and QTL mapping were simultaneously performed using a natural population containing 185 accessions collected worldwide and an F2 population consisting of 240 individuals, respectively. A total of 26 QTLs related to leaf morphological traits were mapped in the F2 population at three different developmental stages, and some QTL intervals were repeatedly detected for different traits and at different developmental stages. Among the 206 significant SNPs identified in the natural population using GWAS, several associated with the leaf thickness phenotype were co-mapped via QTL mapping. By analyzing linkage disequilibrium and transcriptome data from different tissues combined with gene functional annotations, 7 candidate genes from the co-mapped region were identified as the potential causative genes associated with leaf thickness. Conclusions These results presented a valuable cigar tobacco resource showing the genetic diversity regarding its leaf morphological traits at different developmental stages. It also provides valuable information for novel genes and molecular markers that will be useful for further functional verification and for molecular breeding of leaf morphological traits in crops in the future. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-024-05261-8.


Introduction
Leaves are the main photosynthetic organs of plants.Leaf morphology is closely related to plant photosynthetic efficiency, transpiration, and gas exchange, which thereby affect plant biomass.Leaf development involves initiation of leaf primordium, establishment of leaf polarity, and leaf expansion, which is regulated by various factors including developmental signals and environmental conditions [1].
With the development of Next-Generation Sequencing (NGS) technology, studies on Quantitative Trait Locus (QTL) mapping has made great progress, and the QTLs for many important agronomic and economic traits of crops have been mapped.This makes it possible to study the genetic basis of complex quantitative traits, laying the foundation for map-based gene cloning.QTL mapping of some leaf size-related traits in crops have been studied.Jiang et al. (2022) identified five QTL regions regulating leaf size in Alfalfa (Medicago sativa L.) [2]; By combining RIL (Recombinant Inbred Line) and IF 2 (Immortalized F 2 ) populations, Liu et al. (2017) mapped QTLs for maize leaf width [3]; In addition, QTLs for flag leaf length, flag leaf width and flag leaf area in japonica rice were mapped by linkage analysis and Genome-Wide Association Studies (GWAS) using 2 RILs and a natural population [4].Therefore, combining linkage analysis and GWAS to conduct QTL mining is an effective approach for unraveling the genetic basis of leaf morphological characteristics.
Tobacco is widely used as a model plant for fundamental biological research.Nicotiana tabacum L. (common tobacco) is important economic crop and is commonly classified into several categories based on the type and curing process, i.e.Flue-cured tobacco, Sun-cured tobacco, Burley tobacco, Cigar tobacco and Oriental tobacco.A cigar is made up of three parts: wrapper, binder and filler [5].As the outermost part, wrapper tobacco leaves should be flat, wide and thin, with few visible veins; the optimal leaf vein angle is 45 °~70 °.Therefore, cigar tobacco is a good system to unravel the molecular mechanisms regulating leaf development due to its leaf morphological differences from other varieties.In previous studies, multiple QTLs associated with leaf traits of cigar tobacco have been identified on almost all chromosomes using QTL mapping analysis [6].However, there have been no studies using GWAS, not to mention the combination of these two methods for revealing the genetic basis of leaf morphological traits in cigar tobacco.
In this study, germplasm collection of cigar tobacco containing 185 accessions and an F 2 population consisting of 240 individuals derived from two cigar tobaccos were created and used for whole-genome resequencing.Association analysis and QTL mapping were performed using phenotyping data of leaf morphological characters, including leaf flatness and thickness, lateral vein diameter and leaf vein angle, obtained at different developmental stages for three consecutive years in two places.We subsequently combined the two mapping methods and selected the most significant SNPs for linkage disequilibrium (LD) analysis.By considering the leaf-specific expression and gene functional annotations, 7 candidate genes were predicted to be associated with leaf thickness.Our results provide a genetic basis for subsequent gene function verification and for molecular breeding of leaf morphological traits in cigar tobacco as well as other crops in the future.

Plant materials and phenotype measurements
For the natural population, 340 cigar tobacco accessions collected worldwide were first screened phenotypically as well as using SSR markers.A total of 185 representatives (Supplementary Table 1) were selected and used in this study [7,8].All the 185 cigar tobacco germplasms have been deposited in the National Mid-term Genebank for Tobacco which is a collection of tobacco germplasm resource hosted at the Chinese Academy of Agriculture Sciences.The natural population was planted in Shandong (SD), China (120.45°E,36.38°N) and Sichuan (SC), China (104.16°E,31.13°N) in three consecutive years (2019-2021) with three biological replicates and 20 plants per replicate for each accession.Segedinska Ruca, which is characterized by flat and thick leaves, fine lateral veins and large leaf vein angles, was selected as the male parent P 1 for crossing with the female parent Xiawanna which was designated as P 2 and exhibited the opposite phenotype as P 1 , i.e., crinkly and thin leaves, thick lateral veins and small leaf vein angles.The F 2 population was subsequently obtained and consisted of 240 individuals.Together with P 1 , P 2 and F 1 plants, F 2 population were subsequently planted in SD, China, with three biological replicates and 50 individuals for each replicate.
The leaf morphological phenotypes of the natural population, leaf flatness, leaf thickness, lateral vein diameter and leaf vein angle were measured at flowering time.The parents, F 1 and F 2 populations were characterized for leaf traits at the rosette stage, budding time and flowering time.The phenotyping methods used were described previously [9].Briefly speaking, the largest middle leaf was measured using surface roughness meter, angle ruler, vernier caliper and leaf thickness tester (Model No. YH-1) for leaf flatness, leaf vein angle, lateral vein diameter and leaf thickness, respectively.

Whole-genome resequencing, SNP detection and annotation, and genotyping
Young leaves at the rosette stage were sampled and quickly frozen in liquid nitrogen.DNA was extracted by the SLS method and a sequencing library was constructed using an Illumina TruSeq Library Construction Kit.Whole-genome resequencing (WGS) was performed using fastq files following the default parameters, and the sequences were subsequently aligned to the reference genome of N. tabacum L. var ZY300 using BWA.SNP detection was performed using SAMTools (parameters: rmdup; [10]).ANNOVEAR [11] was used to annotate the detected SNPs.

Genome-wide association analysis
The best linear unbiased prediction (BLUP) was calculated using the R package lme4 (CRAN-Package lme4 (r-project.org)).A mixed linear model (MLM) was used to adjust the population structure [12].Moreover, the Bonferroni correction method was used to derive a genome-wide significance threshold.

Construction of the genetic linkage map
To construct a linkage map, we selected sites fixed for alternative alleles between the two parents.JoinMap 4.1 was used to rank the markers of each linkage group, and linkage maps were drawn using the Perl SVG module [13].QTL mapping was performed using the composite interval mapping (CIM) method [14] implemented in WinQTLCart (version 2.5; [15]).The significance threshold was determined using 1000 permutation tests with MapQTL [16].

KASP tag-based genotyping
The 100-200 bp up-and downstream sequences of the SNPs of interest were used for KASP primer design (Supplementary Table 2).KASP-PCR was performed on a 384-well PCR instrument (Eppendorf Mastercycler pro).The KASP-PCR products were read with a BMG POLARstar Omega scanner.The SNP typing results were analysed using KlusterCaller 3.4.1 software (Kincioscience).The aggregated genotypes near the X or Y axis were alleles linked to the FAM or HEX fluorescent tag sequences, respectively.

Both the natural and F 2 populations are suitable for standard genetic analysis
In the natural population, the leaf morphological traits on average in three consecutive years at two different locations exhibited a normal distribution (Fig. 1), with a large degree of dispersion (Table 1).Therefore, the plant The significance of leaf morphological traits between the two parents was also analysed (Table 2).Except for leaf thickness at the rosette stage, there were significant differences in all leaf traits between the two parents at the three developmental stages (P < 0.01), and all the leaf traits investigated showed a certain range of variation (8.61%~37.03%),indicating that the F 2 population is also a qualified genetic population.Moreover, all the leaf traits exhibited a normal distribution at different developmental stages in the F 2 population, which suggested that these traits are suitable quantitative traits and could be used for further genetic analysis (Fig. 2).

Significant SNPs associated with leaf morphological traits were detected via GWAS
A total of 2,119,142 high-quality SNPs were obtained by WGS of 185 accessions in the natural population and were found to be distributed almost evenly on 24 chromosomes, where one SNP appeared in every 440 bp.The numbers of SNPs of different types are shown in Supplementary Table 3. Manhattan plots (BLUP) and QQ plots for each phenotype generated through GWAS analysis are presented in Fig. 3.The significant SNPs shown in the Manhattan peaks are summarized in Supplementary Table 4. Next, we examined the location of these SNPs to determine whether any loci co-mapped in the QTL intervals were identified in the following QTL mapping analysis.

Construction of a high-density genetic linkage map in cigar tobacco
The whole genome of the F 2 population was resequenced to construct a genetic linkage map of cigar tobacco.A total of 1,747,238 SNP markers were obtained, among which 311,687 were of the "aa × bb" type (Supplementary Table 5).After marker screening and filtering, 2,114 of the remaining ones were used for linkage analysis.In the end, the genetic map constructed contained 1682 SNP markers, with a total genetic distance (GD) of 3107.13 cM and an average inner-marker distance of 1.85 cM.Among the 24 linkage groups (lg), the maximum gap between markers was only 30.12 cM.Overall, the density of the genetic map was high (Table 3; Fig. 4A-B).
We then performed collinearity analysis between the genetic map and physical map of ZY300, as shown in Fig. 4C.It was found that each linkage group had high collinearity with the corresponding chromosome, indicating the high quality of the genetic map constructed.

QTL mapping of leaf morphological traits in the F 2 population
A total of 24 QTL sites, for which LODs ranged from 2.81 to 4.85, were detected in the F 2 population and were associated with the four leaf morphological traits during the three developmental periods studied (Table 4).Several QTL intervals were repeatedly detected for different traits and at different periods.
Among them, the interval from 41.9 to 46.3 cM on Hic_ asm_0, associated with leaf thickness at flowering time,

Co-mapping combining the GWAS and QTL mapping results
The physical positions of the significant SNPs and flanking markers of the QTLs were compared using the ZY300 genome.Most importantly, several significant SNPs were found to be located within the QTL intervals (Table 5).
Twelve SNPs associated with leaf thickness co-localized within the QTL between MK163 and MK211.The consistency of the GWAS and QTL mapping results indicated the importance of co-mapped SNPs as potential genetic loci regulating corresponding phenotypes.

Predicted candidate genes associated with leaf thickness in cigar tobacco
To further narrow down the selected QTL region and screen potential candidate genes, a LD plot of the co-mapped region of chromosome Hic_asm_0 was created.A region containing 428 genes with a size of 20 Mb (Hic_asm_0: 178.144-198.139Mb) was found to be linked (Fig. 5A).To validate the regulatory role of the     co-mapped region regarding the leaf thickness phenotype, a haplotype analysis was further performed.The significant difference of leaf thickness between two haplotypes was shown in Fig. 5B.Among the genes detected, 103 genes exhibited missense variance, of which 47 genes were expressed in leaves 2 times higher than in other tissues (Fig. 5C).
Leaf thickness is often associated with cell size, cell number, relative vacuole volume, and cell wall composition [17][18][19].At the molecular level, the process of leaf development involves coordinated regulation among small RNAs, transcription factors and hormones.Based on the functional annotations and relevant literature, 7 genes were predicted to be candidate genes related to the leaf thickness phenotype (Table 6) and were considered very informative for exploring the molecular mechanism controlling leaf thickness in the future.

KASP tag developed for leaf trait genotyping in cigar tobacco
To test whether some SNPs could be applied for genotyping leaf morphological traits in unknown cigar tobacco materials, two significant SNPs detected via GWAS were selected for designing Kompetitive Allele Specific PCR (KASP) markers.The germplasms separated by the genotypes showed significant difference in the phenotypic data (Fig. 6A-B).In the end, the markers LF11 and LVD07 were developed for leaf flatness and lateral vein diameter, respectively (Supplementary Table 6).To verify the KASP markers, a number of germplasms (Supplementary Table Table 6 Predicted candidate genes associated with leaf thickness in cigar tobacco were used for genotyping.The results showed that both the LF11 and LVD07 primers could distinguish the two homozygous genotypes, with one blue signal converged near the X-axis and the other red signal near the Y-axis (Fig. 6C-D).Therefore, both tags were developed successfully and could be used for leaf morphological trait genotyping in cigar tobacco plants in the future.

Discussion
Cigars, as rolled-leaf tobacco products, require cigar wrappers with flat and thin leaves combined with fine lateral veins and large leaf vein angles.The phenotypic differences between wrapper tobaccos and other varieties make it a good model plant for studying the genetic basis of leaf mophological characters.In this study, we used both GWAS and QTL mapping analyses and co-mapped one QTL region to predict causal genes regulating the leaf thickness trait.
Compared with those of other crops, gene mapping studies of tobacco started relatively late, mostly focusing on disease resistance and yield [20,21].Our study is the first to identify the loci and candidate genes associated with the leaf morphological trait in cigar tobaccos.In this study, we preliminarily screened the QTLs responsible for leaf morphological traits through both linkage mapping and GWAS analysis in cigar tobacco.For the leaf thickness phenotype, one QTL within a 30.8Mb interval (Hic_asm_0: 175.087-205.851Mb) from QTL mapping was found co-mapped with the one of 20 Mb (Hic_asm_0: 178.144-198.139Mb) identified by GWAS.Applying of both analysis successfully fine-mapped the QTL intervals, meanwhile increased mapping accuracy.Besides the co-mapped region, some QTLs or SNPs were only detected in one analysis, which suggests either it is a rare allele or no segregation happened at this region in F 2 population.
Fig. 6 The distribution of phenotypic data associated with the markers LF11 (A) and LVD07 (B) and the genotyping results for LF11 (C) and LVD07 (D).The genotypes aggregated near the X or Y axis were alleles linked to FAM or HEX fluorescent tag sequences, respectively Leaf thickness is an important morphological characteristic not only in cigar wrapper breeding but also strongly correlated with plant strategy in response to environmental variables such as high light and CO 2 concentrations [22,23], cold [24][25][26] and drought [27].However, the molecular mechanism regulating leaf thickness is poorly understood.In dicots, leaf thickness is related to mesophyll cell size, cell number and the number of cell layers [28].Thicker leaves tend to have more chloroplasts per cell and more stroma and thylakoid membranes, all of which determine leaf photosynthetic capacity [29].
TARGET OF RAPAMYCIN (TOR), which encodes a highly conserved Ser/Thr kinase, has been reported to regulate leaf cell size as a positive regulator in Arabidopsis [30,31].TOR might directly control cell cycle regulators in leaf development, and surprisingly, cyclin B (CYCB), which partially constitutes the plant cell cycle machinery, was also detected from the co-mapped QTL.The cycb mutant exhibited a reduced leaf area, possibly due to decelerated cell cycle progression [32].In addition to direct influences on cell size and cell number, there are several regulators involved in chloroplast division that ultimately affect cell expansion and cell division.Disruption of FAR1-RELATED SEQUENCE4 (FRS4), a transposase-derived transcription factor, results in enlarged chloroplasts through transcriptional regulation of ACCUMULATION AND REPLICATION OF CHLO-ROPLASTS5 (ARC5), which is essential for chloroplast division [33,34].Furthermore, leaf thickness has been reported to be positively related to cell wall thickness [35].Therefore, several genes related to the metabolism of cell wall components, such as pectinesterase inhibitor 18-like, Extensin-3-like and Cytochrome P450 71A1-like involved in the lignin biosynthetic pathway, were also selected for further functional verification.

Conclusions
In summary, the leaf morphological traits of cigar tobacco, especially leaf thickness, were studied via two combined forward genetic methods, GWAS and QTL mapping.Nine significant SNP peaks and 24 QTL sites associated with four leaf morphological traits were detected in the accessions and F 2 population, respectively.A comparison of the two sets of results revealed a co-localized leaf thickness-related genomic region containing 103 genes which had the missense changes.Through gene expression analysis and interpretation of the gene functional annotation, 7 genes were predicted to be potential functional genes controlling the leaf thickness phenotype during leaf development in tobacco plants.In addition, two KASP markers linked with leaf flatness and lateral vein diameter were developed and verified in some cigar germplasms.These results will lead to the identification of novel genes and molecular markers that will be valuable for further functional verification and for molecular breeding of leaf morphological traits in crops in the future.

Fig. 1
Fig. 1 Distribution of phenotypic data regarding leaf flatness (A), leaf thickness (B), lateral vein diameter (C) leaf vein angle (D) in the natural population had the largest phenotypic contribution (80.25%), indicating its potential as a determining factor regulating the leaf thickness phenotype.One interval, located between 38.4 and 40.1 cM on Hic_asm_0 and associated with leaf thickness, was co-mapped at both the rosette and flowering times.The QTL 13.7 cM ~ 27.4 cM on Hic_asm_12, associated with the leaf vein angle phenotype, overlapped with the QTL 26.4 ~ 34.0 cM on the same chromosome and was linked with the leaf flatness phenotype, indicating that there may be a pleiotropic QTL.

Fig. 2
Fig. 2 Distribution of phenotypic data in the F 2 population.Leaf flatness (A, B, C), leaf thickness (D, E, F), lateral vein diameter (G, H, I) and leaf vein angle (J, K, L) were measured at the rosette stage, budding time and flowering time, respectively

Fig. 4
Fig. 4 High-density genetic map of cigar tobacco (A-B) and collinearity analysis between the genetic and physical maps (C).(A-B) The mapped SNP markers are designated as horizontal blue lines (A) and black lines (B).The X-and Y-axes show the linkage group ID and genetic distance, respectively.(C) The genetic map is designated as red lines, and the physical map is designated as blue

Fig. 5
Fig. 5 LD plot (A) and halpotype analysis (B) of the co-mapped QTL region, and heatmap of leaf thickness-related gene expression in various tobacco tissues (C).(A) Genomic region of chromosome Hic_asm_0:178.144-198.139Mb for leaf thickness is shown in a Manhattan plot and LD heatmap.(B) Haplotype analysis of the co-mapped QTL related to leaf thickness is shown in a violin plot.(C) Expression data were obtained from a publicly available database.Rows represent 103 co-mapped genes that contained missense variance

Table 1
Phenotypic data analysis of leaf morphological traits in the natural population

Table 2
Phenotypic data analysis of leaf traits in the F 2 population Note: The different capital letters represent significant differences between P 1 and P 2 (P < 0.01)

Table 3
The genetic map information

Table 4
The significant QTL intervals detected in the F 2 population

Table 5
The significant SNPs located within the QTL intervals